%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adiabatic sweep through Zeeman structure, Schroedinger & Bloch ansatz.
%
% --- |3>
% --- |2>
% --- |1>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch
    figure(1);				clf;

    % Schrödinger approach approach

    o  =           2*pi*8e3;                                                % Rabi frequency
    t  =   [-1:.001:1]*1e-3;                                                % time
    d  = t/max(t)*2*pi*50e3;                                                % radio frequency sweep
    g  =           2*pi*200;                                                % decay rate
    psi=[1 0 0]';                                                           % initial population

    for kk=1:length(t)
        psi = Hamilton_K(psi,t(2)-t(1),o,o,d(kk),2*d(kk));
        R1(kk)=psi(1)*psi(1)';                                              % populations
        R2(kk)=psi(2)*psi(2)';
        R3(kk)=psi(3)*psi(3)';
    end

%     subplot(223),plot(t*1e3,d,'b','LineWidth',1.2);
%     set(gca,'Linewidth',1.4,'FontSize',16);
%     xlabel('\fontsize{16}{\itt}');
%     ylabel('\fontsize{16}{\it\Delta}   (kHz)');
    subplot(221),plot(t*1e3,R1,'r',t*1e3,R2,'g',t*1e3,R3,'b','LineWidth',1.2);
    set(gca,'Linewidth',1.4,'FontSize',16);
    xlabel('\fontsize{16}{\itt}   (ms)');
    ylabel('\fontsize{16}{\it\rho_{kk}}');
    text(-1.4,1,'\fontsize{16}(a)')

    % Bloch approach
    
    rho=[1 0 0 0 0 0 0 0 0]';                                               % initial population
    for kk=1:length(t)
        rho = Liouville_K(rho,t(2)-t(1),g,g,0,o,o,d(kk),2*d(kk));
        R1(kk)=rho(1);                                                      % populations
        R2(kk)=rho(2);
        R3(kk)=rho(3);
    end

    subplot(222),plot(t*1e3,R1,'r',t*1e3,R2,'g',t*1e3,R3,'b','LineWidth',1.2);
    set(gca,'Linewidth',1.4,'FontSize',16);
    xlabel('\fontsize{16}{\itt}   (ms)');
    ylabel('\fontsize{16}{\it\rho_{kk}}');
    text(-1.4,1,'\fontsize{16}(b)')
end

function psi = Hamilton_K(psi0,dt,o12,o23,d12,d23)
% Schrödinger equation

    H=[
        0         o12/2       0
        o12/2     d12         o23/2
        0		  o23/2       d12+d23
    ];

    psi = expm(1i*H*dt)*psi0;
end

function rho = Liouville_K(rho0,dt,G12,G23,G13,o12,o23,d12,d23)
% Bloch equation

    I=          1i/2;
    g12=         G12;                                                       % decay of coherences
    g13=     G13+G23;
    g23= G12+G13+G23;
    d13=     d12+d23;
    L12=g12/2+1i*d12;
    L13=g13/2+1i*d13;
    L23=g23/2+1i*d23;
    
    M=[
        0		G12    G13    I*o12   -I*o12	0		0		0       0
        0	   -G12    G23   -I*o12	   I*o12	0		0		I*o23	-I*o23
        0		0	  -G13-G23  0		0		0		0		-I*o23	I*o23
        I*o12  -I*o12	0   	-L12	0		I*o23	0		0		0
        -I*o12  I*o12	0       0		-L12'	0		-I*o23	0		0
        0		0		0		I*o23	0		-L13	0		-I*o12	0
        0		0		0		0		-I*o23	0		-L13'	0		I*o12
        0       I*o23  -I*o23	0		0		-I*o12	0		-L23	0
        0      -I*o23	I*o23	0		0		0		I*o12	0		-L23'
    ];

    rho = expm(M*dt)*rho0;
end
